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FOREWORD 

This  paper  was  prepared  for  presentation  at  the  Air  Force  Special 
Weapons  Center  Symposium,  ''Instrumentation  for  Nuclear  Weapons  Effects 
Simulation,"  12-13  March  1970. 

The  filter  set  described  in  this  paper  was  developed  during  normal 
operation  of  the  Physical  Sciences  Branch,  Nuclear  Weapons  Effects  Divi¬ 
sion,  U.  S.  Army  Engineer  Waterways  Experiment  Station  (WES),  under  the 
general  supervision  of  Mr.  G.  L.  Arbuthnot,  Jr.,  Chief,  Nuclear  Weapons 
Effects  Division,  and  Mr.  L.  F.  Ingram,  Chief,,  Physical  Sciences  Branch. 

Mr.  II.  D.  Carleton  was  responsible  for  project  development  and  documenta¬ 
tion.  The  filters  of  the  set  were  incorporated  into  computer  program 
number  803-G9R0-118  during  the  period  from  June  to  December  1969  by 
Mr.  J.  T.  Brogan,  assisted  by  Mrs.  C.  ,J.  Pittman. 

This  work  was  performed  to  improve  the  quality  of  ground  shock  data 
obtained  on  field  tests  sponsored  largely  by  the  Defense  Atomic  Support 
Agency. (DASA) .  Additional  funding  support  was  provided  from  DASA  Project 
61102H-L11CAXSX502,  "Ground  Stress  and  Motion  Measurements."  Computer  pro¬ 
grams  developed  as  a  result  of  this  work  are  also  being  used  for  the  proc¬ 
essing  of  laboratory  test  data  obtained  from  various  nuclear  weapons 
effects  research  (NWER)  studies  in  the  WES  blast  load  generator  facility. 

WES  Director  during  the  filter  set's  development  was  COL  Levi  A. 
Brown.  Technical  Director  was  Mr.  F.  R.  Brown. 
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NOTATION 


til 

a  the  u  n  coefficient  in  a  filter  numerator  polynomial 

U  til 

b  the  v  coefficient  in  a  filter  denominator  polynomial 

/ 

e  natural  logarithm  base 

f  any  desired  frequency  (hertz) 

f^  any  aliased  frequency  (hertz) 

f^  folding  or  Nyquist  frequency  (hertz) 

f  sampling  frequency  (hertz) 

s 

j  the  imaginary  operator  V^T 

K  maximum  power  of  z  in  the  numerator  of  a  filter  function 

L  maximum  power  of  z  in  the  denominator  of  a  filter  function 

n  sample  number  (an  integer) 

s  the  generalized  frequency  variable  of  Laplace 
t  elapsed  time  (seconds  or  sampling  units.) 

T  sampling  period  (seconds) 

x  wave  displacement  value  for  any  considered  input  time  =  input 
sample  number  n 

xn+1  the  input  sample  immediately  following  (in  time)  input  sample 
number  n 

xn  1  the  input  sample  immediately  preceding  input  sample  number  n 

xn  2  the  input  sample  immediately  preceding  input  sample  number  n  -  1 

yn  the  output  sample  corresponding  to  input  sample  number  n 
y^  ^  the  output  sample  immediately  preceding  output  sample  nymber  n 

yn  2  the  output  sample  immediately  preceding  output  sample  number  n  -  1 

z  the  z  transform  variable 

jzj  the  modulus  of  z 

Af  any  difference  in  frequency  (hertz) 

vii 


T]  the  real  component  of  z 
v  the  imaginary  component  of  z 
a  the  real  component  cf  s 

u)  the  frequency  variable  (radians  per  unit  time),  and  imaginary 
component  of  s 
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SUMMARY 


Digital  tapes  produced  from  analog  field  data  often  carry  over  severe 
noise  problems.  A  batch  data  processing  capability  should  therefore  be  com¬ 
plemented  by  digital  filters  which  can  remove  noise  routinely.  This  paper 
describes  a  set  of  filters  developed  for  use  with  Nuclear  Weapons  Effects 
Division  standard  data  processing  codes.  Emphasis  in  assembly  of  the  set 
was  placed  on  accuracy  of  signal  retention  and  on  adaptability  to  general 
purpose  computers. 

Part  I  discusses  data  which  is  sampled  at  equally  spaced  times,  the 
use  of  sample  numbers  to  designate  time  positions,  and  the  requirement  for 
use  of*  a  sampling  frequency  at  least  eight  times  greater  than  the  highest 
expected  signal  frequency. 

Parts  II,  III,  and  IV  document  the  filters  of  the  three  subroutines 
presently  being  used  for  noise  removal.  The  first  subroutine  eliminates 
randomly  spaced  single-sample  ''spikes''  through  the  application  of  inequality 
conditions  to  data.  The  second  subroutine  provides  frequency  filters  for 
the  removal  of  zero  drift  and  single  noise  frequencies.  Low- pass  frequency 
filters  are  available  in  the  third  subroutine. 

Part  V  discusses  the  ten  example  plates  used  to  demonstrate  the  filter 
set’s  effectiveness  and  speed  of  operation. 


DIGITAL  FILTERS  FOR  ROUTINE  DATA  REDUCTION 


PART  I:  DIGITAL  DATA  SAMPLING 


Introduction 

1.  With  the  purchase  of  an  off-line  analog  to  digital  (a/d)  converter 
in  1968,  the  Nuclear  Weapons  Effects  Division  (NWED)  of  the  U.  S.  Army 
Engineer  Waterways  Experiment  Station  (WES)'  obtained  a  batch  data  process¬ 
ing  capability.  Since  analog  data  is  often  quite  noisy,  and  because  the 
A/D  converter  changes  both  signal  and  noise  to  digital  form,  digital  tech¬ 
niques  for  noise  removal  became  an  immediate  consideration.  Studies  of 
several  approaches  to  the  noise  problem  led  to  the  conclusion  that  a  digi¬ 
tal  filter  set  should  be  assembled  for  routine  use.  Filters  accepted  for 
the  set  were  to  cause  no  phase  distortion  among  passed  frequencies  and  to 
operate  with  reasonable  econony  on  the  general  purpose  computer  equipment 
available  at  WES. 

2.  This  report  has  been  written  to  document  the  filter  set  developed, 
and  to  aid  project  officers  in  their  use  of  its  component  filters.  These 
filters  are  included  as  subroutines  to  NWED  standard  data  processing  codes. 

A  basic  familiarity  with  complex  numbers  is  assumed  for  the  discussions 
contained  in  Parts  III  and  IV. 

Notes  on  Discrete  Sampling 

3.  If  the  displacement  of  a  continuous  wave  form  is  sampled  periodi¬ 
cally,  the  resulting  time,  sequence  of  equally  spaced  observations  is  said 
to  be  a  "discrete  time  series."  The  time  represented  by  any  given  sample 
in  this  series  would  be  t  =  nT  ,  where  n  is  the  sample  number  (counting 
by  units  from  sample  number  0  at  t  =  0)  and  T  is  the  sampling  period  (a 
constant ,.  commonly  in  seconds,  which  is  the  reciprocal  of  the  sampling 
frequency  fg  ,  commonly  in  hertz) .  If  T  is  defined  as  one  unit  of  time, 
t  =  nT  becomes  t  =  n  . 
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4.  Fig.  1  shows  examples  of  periodically  sampled  continuous  waves. 

At  a  given  sampling  frequency,  say  for  example  6000  Hz,  no  oscillation 
can  be  represented  if  all  sampled  displacement  values  are  exactly  the  same 
(as  in  fig.  la) .  This  is  because  actual  displacement  between  samples 
cannot  be  determined  by  examination  of  the  discretely  sampled  data.  Only 
estimation  or  interpolation  is  possible.  If  ail  displacement  sample  values 
are  exactly  the  same,  and  not  equal  to  zero,  "zero  frequency"  is  repre¬ 
sented  by  the  constant  offset  from  the  zero  displacement  base  line  (dis¬ 
placement  equilibrium  position) .*  If  only  two  displacement  values  exist 
among  all  the  samples  taken,  and  these  two  values  are  found  alternating 

at  every  change  in  sample  through  time  (as  in  fig.  lb),  then  the  frequency 
represented  by  this  alternation  is  half  the  sampling  frequency,  or  for  our 
example,  3000  Hz.  This  frequency,  known  as  the  folding  or  Nyquist  fre¬ 
quency  fN  ,  is  the  highest  frequency  which  may  be  represented  by  discretely 
sampled  data  at  a  given  sampling  frequency.  Frequencies  higher  than  fold¬ 
ing,  if  they  are  present  on  the  original  continuous  data,  will  be  repre¬ 
sented  by  distortion  of  frequencies  below  folding..  Any  frequency  f  =  + 

Af  ,  where  0  <  Af  <  f N  ,  will  be  "seen"  by  the  sampling  process  as  the 
"folded"  or  "aliased"  frequency,  f  A  =  f  N  -  Af  .  For  the  example  sample  rate 
of  6000  Hz,  for  instance,  4200-Hz  noise,  if  present,  would  be  "frequency 
aliased"  as  1800  Hz  and  added  to  any  existing  l800-Hz  signal,  since  folding 
frequency  would  be  3000  Hz.  (Additional  discussion  of  aliasing  may  be  found 
in  Blackman  and  Tukey,-'-  pp  216-219  and  521-524.) 

5.  It  has  been  stated  in  the  preceding  paragraph  that  folding  fre¬ 
quency  is  the  highest  frequency  which  may  be  represented  by  discretely 
sampled  data  at  a  given  sampling  frequency.  A  reasonably  accurate  approxi¬ 
mation  of  folding  frequency  sine  waves  can  be  made  only  where  the  times  of 
sampling  correspond  to  the  times  of  maximum  displacement  for  a  wave  (as  in 

s 

fig.  lb).  If  the  times  of  maximum  wave  displacement  are  out  of  phase  with 
the  sampling  times,  both  displacement  and  phase  angle  are  unacceptably  al¬ 
tered  for  the  wave  represented  (see  fig.  2a) .  In  fact,  a  90-deg  phase  dif¬ 
ference  between  sample  times  and  maximum  displacement  times  for  the  folding 

*  Constants  and  linear  trends  are  zero  frequency.  See  Blackman  and 
Tukey,-^-  section  19. 
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Fig.  1.  Periodically  sampled  continuous  waves 
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Fig.  2.  Distortion  of  sine  waves  due  to  sampling  and  plotting 


i 


frequency  reduces  all  sample  displacements  to  zero,  and  no  wave  is  repre¬ 
sented  (see  fig.  2b).  Frequencies  below,  but  near,  folding  are  also 
strongly  affected  by  the  phase  relation  between  sampling  and  maximum  dis¬ 
placement  times.  The  amplitude  response  which  may  be  expected  over  por¬ 
tions  of  a  wave  with  sampling  times  at  the  least  favorable  positions  may  be 
determined  by  finding  the  cosine  of  the  angle  (l80  fT)°,  where  fT  is  the 
reciprocal  of  the  number  of  samples  per  cycle  for  the  frequency  f  being 
considered  (see  figs  2b,  c,  and  d) .  Fig.  3  is  a  plot  of  amplitude  responses 


for  sampled  frequencies  at  the  most  unfavorable  sampling  phase  relations 


Sampling  Frequency  Selection 

6.  The  selection  of  a  given  sampling  frequency  band  limits  signal 
frequency  representations  by  establishing  folding  frequency  at  half  the 
frequency  of  sampling.  The  plotting  process,  through  its  inability  to 
superimpose  a  sine  wave  upon  all  displacement  samples  which  may  have  been 
taken  from  higher  frequency  sine  wayes,  further  limits  the  frequencies 

V. 

which  may  be  accurately  represented  upon  output  plots  (as  shown  in  figs.  2 
and  3).  Sampling  frequency  should  be  at  least  eight  times  greater  than 
the  highest  expected  signal  frequency,  to  assure  reasonable  plot  accuracy. 
Data  sampled  at  this  minimum  rate  will  be  relatively  economical  to  proc¬ 
ess.  Higher  sampling  frequencies  may  be  necessary  where  the  highest 
signal  band  frequencies  must  be  reproduced  with  extreme  accuracy,  or 
where  noise  frequencies  will  be  aliased  into  the  signal  band  by  lower 
sampling  frequencies . 


PART  II:  THE  RANDOM-SPIKE  REJECTION  FILTER 


7.  Digital  tapes  occasionally  have  upon  them  single  samples  (at  ap¬ 
parently  random  locations  in  the  discrete  time  series)  which  show  very  high 
displacement  values,  and  which  hear  no  true  relation  to  the  seismic  data 
thereon.  These  spurious  "spikes"  frequently  cause  large  errors  during  fre¬ 
quency  filtering,  integration,  and  frequency  analysis  operations.  The 
problem  may  be  eliminated  by  application  of  mathematical  inequalities  to 
the  raw  digital  data. 

8.  The  main  condition  which  defines  a  spike  is: 


1.25 


x  . 

>  o-75 

n+X 


where  xn  represents  the  value  of  the  input  (displacement)  presently 

under  consideration,  x  ..  represents  the  immediately  preceding  input 

n-i 

value,  and  x^+1  represents  the  next  input  value  in  the  future.  The  above 
inequality  simply  states  that  the  ratio  of  the  two  "sides"  of  a  spike  must 
be  nearly  unity. 

9.  The  condition  which  defines  the  sizes  of  spikes  to  be  removed  is: 


,  ■  ow*) 


>  220 


In  this  case  spikes  larger  than  220  units  are  defined  as  spurious. 

10.  Instructions  to  the  computer  tell  it  to  consider  each  sample  in  a 
discrete  time  series  against  the  above  inequalities.  If  any  sample  (rela¬ 
tive  to  the  samples  on  either  side  of  it)  meets  both  of  the  above  condi¬ 
tions,  it  is  a  spurious  spike,  and  the  computer  is  instructed  to  replace  it 

/x  i  +  xn+1  \ 

with  a  sample  displacement  value  of  l — - — ^ - 1 ,  which  is  simply  the 


average  of  the  adjacent  values.  The  only  samples  which  will  be  affected  in 
any  way  by  this  program  are  those  single  sample  spikes  larger  in  size  than 
the  specified  minimum.  This  minimum,  of  course,  may  be  carefully  altered 
where  it  appears  from  outputs  that  the  program  is  not  operating  at  best 
efficiency. 
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PART  III:  SELECTIVE  FREQUENCY  REJECTION 


Complex  Plane  Frequency  Representations 


11,  The  discussion  in  this  section  relates  graphical  frequency  rep¬ 
resentations  to  Fourier  series  theory  and  Laplace's  generalization  of  this 
theory.  (See  a  text  such  as  LePage  for  background.  Additional  informa¬ 
tion  may  he  found  in  Blackman  and  Tukey,  pp  252-258.) 

12.  The  Fourier  frequency  variable  u>  (expressed  in  radians  per 

unit  time)  is  related  to  Laplace's  s  ,  a  generalization  of  ou  ,  as  fol- 
sT  i(i)T 

lows :  e  is  a  generalization  of  ed  (where  T  is  the  sampling 
period,  a  constant)-,  and  s  =  a  +  jcu  ,  where: 

s  is  called  the  "generalized  frequency  variable" 

a  is  the  real  component  of  s 

j  is  the  imaginary  operator 

to  is  the  frequency  variable,  and  the  imaginary 
component  of  s 

Therefore,  s  is  a  complex  number,  the  value  of  which  is  represented  by  the 

ordered  pair  of  real  numbers  (a,uu)  (see  fig.  4).  While  itself  real,  u>  is 

used  in  this  terminology  as  the  multiplier  of  j  in  the  s  plane.  When 

sT  *ii)T 

ct  =  0  ,  s  =  Ju)  ,  and  e  =  e^  =  z  .  In  this  specific  case,  values  cf 


Fig.  4.  s  plane  representation  of  frequency 
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z  vary  as  cu  changes,  since  T  is  a  constant. 

13.  The  s  plane  is  used  to  define  values  of  s  ,  the  generalized 

sT 

frequency  variable.  If  the  function  e  of  the  independent  complex 

variable  s  is  plotted  on  a  second  complex  plane  (the  z  plane)  as  s 

follows  a  path  along  the  s  plane  imaginary  axis  from  (0,-rr)  to  (0,tt), 

the  unit  circle  (modulus,  jzj  =  l)  will  result  (see  fig.  5)*  This  map- 

iooT 

ping  is  done  through  the  use  of  Euler’s  identity,  e°  =  cos  uff  +  j  sin  u/T  , 
where  u)T  is  the  phase  angle  on  the  z  plane.  This  angle  represents 


Fig.  5*  z  plane  representation  of  frequency 

that  portion  of  a  full  cycle  or  period  at  any  frequency  f  which  would 
be  completed  during  one  sampling  period  T  .  Half  the  unit  circle  (n 
radians),  then,  represents  folding  frequency  f^  since  folding  is 
half  of  the  sampling  frequency,  and  only  half  of  a  cycle  would  be 
completed  during  one  sampling  period.  In  terms  of  folding  frequency: 

U)T  =  2iTfT  =  2rrf  ~  (l) 

N  N 

Zero  frequency  is  represented  by  z  =  (1,0)  ,  and  folding  frequency  by 
z  =  (-1,0)  .  Values  of  frequency  between  zero  and  folding  are  represented 
by  u>T  intercepts  between  0  and  tt  on  both  the  upper  and  lower  halves  of 
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the  unit  circle,  reflecting  the  band  limited  nature  of  frequency  spectra 
for  discrete  time  series.  Complex  conjugate  values  of  z  (conjugate  val¬ 
ues  frcm  above  and  below  the  real  axis)  are  used  to  produce  real  coeffi¬ 
cients  in  the  synthesis  of  filter  polynomials.  If  T  is  taken  as  unity, 
a)  becomes  equivalent  to  the  phase  angle  on  the  z  plane.  (Additional 
discussion  of  frequency  representations  may  be  found  in  Robinson  and 

O 

Treitel.  ) 


Filter  Synthesis 

14.  An  important  usefulness  of  frequency  representation  in  the  z 
plane  lies  in  its  adaptability  to  filter  synthesis i  Digital  filters  may  be 
expressed  as  a  ratio  of  two  polynomials  in  z  ,  where  the  roots  of  the 
numerator  are  zeros,  and  those  of  the  denominator,  poles  (Shanks,^  pp  35- 
4l).  For  simple  single  frequency  rejection  filters,  zeros  and  poles  may 
be  determined  by  examination  of  the  unit  circle  in  the  z  plane .  Because 
frequency  rejection  is  desired,  values  of  u)  (represented  by  z  plane 
points  on  the  unit  circle)  must  cause  the  filter  function  (the  ratio  of  two 
polynomials)  to  go  to  zero  at  the  appropriate  frequency.  This  is  done  by 
the  choice  of  a  zero  or  zeros  on  the  unit  circle.  All  other  values  of  u>  , 
however,  must  produce  a  filter  function  response  as  near  as  possible  to 
unity,  since  distortion  of  nonrejected  frequencies  is  not  desirable.  The 
selection  of  a  pole  or  poles  just  outside  of  the  mit  circle  at  the  same 
phase  angle(s)  as  the  zero(s)  will  produce  this  result,  at  the  same  time 
keeping  the  function  stable.  (Because  z  plane  representations  of  m 
(where  a  =  0)  are  confined  to  the  unit  circle,  the  filter  function  will 
not  become  infinite  at  any  considered  point.) 

The  Zero-Frequency  Rejection  Filter 

15.  An  offset  oscillating  time  function  can  be  made  to  seek  the  zero 
displacement  base  line  as  its  equilibrium  position  by  application  of  a 
zero-frequency  digital  filter  designed  by  the  method  outlined  in  the  pre¬ 
ceding  paragraph.  As  shown  previously,  zero  frequency  is  represented  by 
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the  point  z  =  (1,6)  on  the  unit  circle.  The  linear  polynomial  having 
1  +  jO  as  its  root  is  1  -  z  ,  which  will  he  the  numerator  of  the  zero- 
frequency  rejection  filter  function.  Choosing  a  pole  value  on  the  cu  =  0 
radial  hut  just  outside  of  the  unit  circle,  we  have  z  =  (1.01,0).  The 
linear  polynomial  having  1.01  +j0  as  its  root  is  1.01  -  z  ,  the  denomi¬ 
nator  of  the  zero-frequency  reject  filter  function.  The  filter  function, 
then,  is: 


■  raBr  <2> 

1 6.  A  recursive  algorithm  (Shanks,^  pp  3^-35)  may  he  used  to  apply 
this  filter  function  to  an  input  time  series.  The  general  recursion  equa¬ 
tion  for  rational  filters  is: 

K  L 

yn  =  ^  Vn-u  -  2  Vn-v  (3) 

u=0  v=l 

Where  a  values  represent  the  coefficients  in  the  numerator  polynomial, 
h  values  represent  the  coefficients  in  the  denominator  polynomial, 
x  values  represent  input  displacement  samples,  and  y  values  represent 
output  displacement  samples.  The  values  of  u  and  v  are  taken  from  the 
powers  of  z  associated  with  the  coefficients  of  the  filter  function  poly¬ 
nomials,  and  n  is  the  sample  number  of  the  output  displacement  sample 
being  computed.  K  is  the  maximum  power  of  z  represented  in  the  numera¬ 
tor,  and  L  the  maximum  power  of  z  found  in  the  denominator.  The  coef¬ 
ficient  of  the  first  term  in  the  denominator  hQ  should  he  unity.  Where 
this  is  not  the  case  for  a  particular  filter  function,  both  numerator  and 
denominator  are  divided  by  the  initial  value  of  hQ  (excepting  the  case 

where  h  =  0)  to  bring  b  to  unity, 
o  o 

17.  In  the  case  of  the  zero-frequency  rejection  filter  function, 
numerator  and  denominator  must  be  divided  by  1.01,  which  gives: 


F(z)  =  0; -  0.??00??z 
'  '  1  -  0.990099z 


(4) 


By  application  of  the  recursion  summation  (equation  3) : 
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yn  =  0.990099xn  -  0.990099*^  +  0.990099yn.1 

yn  =  0.990099(xn  -  xn_1  +  yn_x)  (5) 

Equation  5  may  be  used  directly  in  an  appropriate  computer  program  to  re¬ 
move  the  zero-frequency  component  from  a  discrete  time  series.  The  recur¬ 
sion  method  of  application  for  this  filter  has  the  advantages  of  speed  and 
accuracy  as  compared  to  other  methods  of  filter  application,  because  the 
consideration  of  a  previously  computed  output  '(y  and  one  previous  in¬ 
put  value  (xn-1)  for  each  output  computation  (y  )  eliminates  the  need  for 
consideration  of  a  very  long  (theoretically  infinite)  series  of  previous 
inputs  for  each  computed  output  value.  A  single  operation  (called  a 
"pass")  with  this  filter  brings  an  offset  wave  to  the  base  line  exponen¬ 
tially.  Thus,  early  displacement  values  in  the  time  series  are  only 
slightly  shifted,  and  the  once-filtered  wave  may  appear  to  have  a  "warp" 
similar  to  that  of  a  very  low  frequency  harmonic.  A  second  operation,  this 
time  on  the  time  reversed  output  of  the  first  pass,  will  correct  this 
situation  (Shanks,^  pp  41-42).  The  phase  response  of  this  "two-pass"  fil¬ 
ter  will  be  zero  for  all  frequencies.  The  amplitude  response  is  shown  in 
fig.  6. 


The  Single -Frequency  Rejection  Filter 


18.  Where  a  single  frequency  has  been  identified  as  spurious,  it  may 
be  removed  from  a  discrete  time  series  by  an  operation  only  slightly  more 
involved  than  that  used  for  the  zero-frequency  rejection  filter.  From 
equation  1  (with  the  sampling  period  T  taken  as  unity),  it  is  seen  that 
the  angle  tu  which  represents  any  temporal  frequency  f  is  dependent  on 
the  ratio  of  that  frequency  to  the  folding  frequency.  Recalling  that  fold¬ 
ing  frequency  is  equal  to  half  of  the  sampling  frequency: 


U)  = 


2rr f 
f 

s 


(6) 
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Say  that  the  frequency  to  be  rejected  is  960  Hz,*  and  that  sampling  fre¬ 
quency  is  12,000  Hz. 

u>  =  =  0.5026548  radians  =  28.8  deg 

The  points  of  intercept  of  u)  (plus  and  minus)  with  the  unit  circle  are 
zeros.  The  projections  of  these  points  on  the  two  axes  of  the  complex 
z  plane  (fig.  5)  are: 

T]  =  +1  cos  u)  =  +cos  28.8  deg  =  +O.8763067 


v  =  +1  sin  a)  =  +sin  28.8  deg  =  +0.4817537 
For  p)les  | z  |  =  1.01  on  the  same  radials  as  the  zeros: 

71  =  +1.01  cos  u,  =  +1.01  cos  28.8  deg  =  +0.8850697 


v  =  +1.01  sin  u)  =  +1.01  sin  28.8  deg  =  +0.4865712 

Therefore,  the  zeros  and  poles  of  a  digital  filter  which  removes  960  Hz 
from  data  sampled  at  12,000  Hz  are: 


two  zeros 


two  poles 


/  z  =  0.8763067  +  j  0.4817537 
\  z  =  0.8763067  -  j  0.4817537 

(  z  =  0.8850697  +  j  0.4865712 
(  z  =  0.8850697  -  j  0.4865712 


To  determine  the  filter  function  polynomials,  form  equations  of  the  lowest 
possible  degree  (with  real  coefficients)  which  have  the  above  values  of  z 
for  roots.  The  numerator  polynomial  is  formed  from  the  two  zeros;  the 
denominator  polynomial  from  the  two  poles.  For  this  filter  then: 


*  A  common  noise  frequency— the  result  of  a  60-Hz  source  and  a  16:1 
digitizing  tape  speed  ratio,  as  in  plate  9* 
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F(z)  =  1  -  1»7526134z  +  z2 

1.0201  -  1.7701395z  +  z2 

Dividing  both  numerator  and  denominator  by  1.0201  to  bring  b  =  1  : 

F(z)  -  Q»9802960  -  l.,7l80800z  +  0.980296022 
1  -  1.7352608z  +  0.9802960z2 

Putting  F(z)  in  recursive  form  by  the  use  of  equation  3: 

yn  =  0.9802960xn  -  1.7180800x  x  +  0.9802960xn_2 
+  1.7352608yn_1  -  0.9802960yn_2 

Because  the  above  filter  function  polynomials  are  quadratic,  the  recursive 
form  of  this  filter  must  consider  input  and  output  samples  up  to  two  sample 
periods  previous  to  the  present  input.  It  nevertheless-  retains  great  speed 
of  computation  by  comparison  with  other  methods  of  filter  application. 
Applied  by  the  two-pass  method,  it  is  a  zero  phase  filter.  Amplitude  re¬ 
sponse  is  shown  in  fig.  7. 

19.  It  is  not  necessary  to  go  through  all  of.  the  foregoing  steps  in 
the  synthesis  of  a  single -frequency  reject  filter.  The  recursive  equation 
for  any  reject  frequency/sampling  frequency  combination  for  tu  will  always 
take  the  form: 


yn  =  0.9802960xn  -  Axn_1  +  0. 98029 60xn_2  +  By^  -  0.9802960yn_2 


where 


O  /2Trr  \ 

2  008  (-) 
1.0201 


B  = 


2.02  cos 


1.0201 


As  an  evidence  of  the  exponential  nature  of  this  algorithm,  recall  from 
Euler's  identity  that  it  can  be  shown  that: 


cos  gu 


2 


f 

Y~  =  frequency  as  a  fraction  of  sampling  frequency 
s 

.  7.  Amplitude  response  of  single-frequency  rejection  filter 


PART  IV:  LOW- PASS  FILTERS 


Filter  Synthesis 


20.  Low-pass  filters  are  a  virtual  necessity  where  many  high  noise 

frequencies  accompany  shock  data.  Obviously,  a  filter  function  must  have 

numerous  zeros  if  the  represented  filter's  response  is  to  be  kept  very  near 

zero  over  a  broad  band  of  frequencies.  This  requirement  can  lead  to  filter 

5 

functions  with  many  terms  and  coefficients .  Whittlesey  suggests  an  alter¬ 
native  to  the  use  of  multicoefficient  filter  functions.  Simple  filter  func¬ 
tions  may  be  used  in  series  to  form  composite  low-pass  filters.  This  ap¬ 
proach  limits  the  choice  of  low-pass  bands  to  certain  specific  fractions 
of  the  sampling  frequency,  but  filter  operation  becomes  relatively  economi¬ 
cal,  and  the  routine  can  be  two-passed  to  eliminate  phase  distortion. 

21.  Four  composite  low-pass^  filters  may  be  constructed  by  a  proper 
combination  of  the  following  components: 


Operator  Operation  to  be  Performed  to 
Designation  Yield  Each  Displacement  Output 


= 

X 

4* 

X  . 

n 

n-1 

X 

4- 

X  n 

n 

n-2 

X 

JL 

x  » 

n 

n-3 

=: 

X 

4* 

n 

n-4 

— 

X 

4* 

x  c 

n 

n-6 

X 

4* 

X  Q 

n 

n-8 

ST 

X 

+ 

X  , 

n 

n-1 

+  x.  „ 

n-d 


y  =  0.43  (x  -  y  «)  +  x 
Jn  v  n  Jn-3  n 


Filter  Function 
(z  Plane) 

F(z)  =  1  +  z 

F(z)  =  1  +  z2 

F(z)  =  1  +  z3 

F(z)  =  1  +  z^ 

F(z)  =  1  +  z 

F(z)  =  1  -i-  z^ 

r 

F(z)  =  1  +  z  +  z 


F(z)  = 


=  0.36  (x„  -  yA  +  x„  F(z)  = 


y  =  0.32  (x  -  y  n)  +  x  F(z)  = 
Jn  '  n  ^n-9  n  '  ' 

y  =  0.324  (x  -  „•  , J  +  x  F(z)  = 
°n  '  n  "a- 12'  n  '  ' 


1.43 

1  +  0.43z3 
1.36 

1  +  0.36z6 
1.32 

1  +  0.32z9 
1.324 

1  +  0.324z] 


An  operator  with  an  "S"  designation  is  a  summing  operator;  one  with  an 
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22.  The  zeros  of  a  composite  filter  are  the  zeros  of  the  summing 
operators  which  are  included  in  the  makeup  of  the  composite  filter.  A 
recursion  operator  is  used  in  each  composite  filter  to  "square  up"  the  low- 
pass  band  of  that  filter’s  response  curve.  One  gain  (amplitude)  correc¬ 
tion  is  made  to  each  composite  filter’s  output  following  the  series  of 
summing  and  recursion  operations. 

The  Low-Pass  Options 

23.  The  option  1  low-pass  filter  uses  four  operators  in  the  following 
sequence:  SI,  SX,  S2,  R3.  The  operator  SI  is  applied  at  each  displacement 
sample  in  the  wave  to  be  filtered,  to  produce  an  output  sample.  The 
series  of  SI  output  displacement  samples  is  then  treated  as  input  to  the 

SX  operator.  Each  succeeding  output  is  treated  as  input  to  a  following 
operator  until  the  series  is  complete.  An  output  from  the  R3  operator  for 
the  option  1  low-pass  filter  has  been  amplified  twelve  times  by  the  series 
of  operations.  For  this  reason,  each  final  displacement  sample  is  multi¬ 
plied  by  the  gain  factor,  l/l2.  A  second  operation  of  the  entire  filter 
series,  this  time  on  the  time  reversed  output  of  the  first  pass,  brings 
phase  response  for  all  frequencies  to  zero.  The  amplitude  response  for  the 
two-passed  option  1  filter  is  shown  in  fig.  10. 

24.  Each  of  the  remaining  options  uses  operators  and  a  gain  factor 
as  shown  in  proper,  sequence  below: 

Option  2:  SI,  SX,  S2,  S2,  S4,  R 6,  l/48 

Option  3:.  SI,  SX,  S3,  S2,  s6,  R9,  1/48 

Option  4:  SI,  SX,  S2,  S4,  S2,  S8,  R12,  1/96 

Each  filter  is  two-passed  to  eliminate  phase  distortion.  Amplitude  re¬ 

sponses  for  options  2,  3,  and  4  are  shown  in  figs.  11,  12,  and  13, 
respectively. 
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Fig.  11.  Amplitude  response  of  option  2  low-pass  filter 
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PART  V:  USE  OF  THE  FILTER  SET 


General 


25.  The  filters  described,  in  this  paper  depend  for  operation  upon 
data  inputs  which  are  sampled  at  equally  spaced  times.  An  examination  o’f 
the  various  response  curves  shows  that  the  filters*  frequency  characteris¬ 
tics  are  dependent  upon  sampling  frequency.  If  sampling  frequency  is  rela¬ 
tively  low,  response  curve  roll-offs  and  rejection  notches  affect  fewer 
frequencies  than  would  be  affected  if  sampling  were  faster.  Where  signal 
and  noise  frequencies  are  in  close  proximity,  careful  comparisons  of  proc¬ 
essing  requirements  and  response  curves  should  be  made  before  sampling  fre¬ 
quency  is  selected. 

26.  Plates  1-10  show  examples  of  digital  data  improvement  using  the 
filters  described  in  this  paper.  Each  plate  has  been  labeled  to  show  gage 
type,  sampling  frequency,  filters  applied,  and  computer  time  used  in  the 
filters.  Where  integrations  are  included  with  gage  data,  they  are  produced 
from  raw  or  filtered  gage  outputs.  No  filter  has  been  applied  directly  to 
any  integration  output. 


Random-Spike  Rejection 

27. -  Large  single  sample  spikes  affect  data  plots  as  shown  in 
plate  1.  The  computer-set  scales  for  the  accelerometer  plot  create  a  low- 
amplitude  data  presentation  because  of  the  spurious  spikes’  high  amplitude. 
A  second  problem  is  seen  in  the  integration  to  velocity;  the  spikes  have 
distorted  this  curve.  Plate  2  shows  the  same  data  after  operation  of  the 
random  spike  rejection  filter.  Where  used,  this  subroutine  precedes 
frequency  filtering  operations. 

Low-Pass  Operations 

28.  Plates  2,  5,  and  7  show  data  affected  by  high  noise  frequencies, 
while  plates  3>  8,  and  8  show  the  same  data  after  application  of 
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appropriate  low-pass  filters .  The  strain  gage  data  in  plates  7  and  8  were 
sampled  very  rapidly  because  it  was  included  on  a  tape  with  accelerometer 
data.  Use  of  the  option  4  low-pass  filter  was  necessary  in  this  case. 

All  four  low-pass  options  are  included  in  one  NVJED  subroutine. 

Selective  Frequency  Rejection 

29.  The  integration  to  velocity  in  plate  3  clearly  indicates  accumu¬ 
lation  of  area  between  the  accelerometer  plot  and  its  zero  ordinate— i.e. , 
an  offset.  This  offset,  though  very  small,  distorts  the  integration  plots 
which  are  intended  to  represent  particle  velocity  and  displacement.  Plate  4 
shows  the  same  data  after  application  of  the  zero-frequency  rejection  filter 
to  the  accelerometer  output.  This  filter  should  be  used  with  care,  since 
surface  waves,  if  present,  are  likely  to  fall  in  the  rejection  notch.  Also, 
success  with  zero- frequency  rejection  depends  on  signal  oscillation  which 
centers  approximately  on  its  trace’s  equilibrium  position.  Large  offsets 
are  handled  by  making  an  estimated  base-line  shift  before  filtering. 

30.  The  960-Hz  periodic  wave  on  the  stress  signal  in  plate  9  has  been 
removed  in  plate  10  by  application  of  the  single-frequency  rejection  filter. 
Note  retention  of  the  sharp  peak  at  1.2  msec.  This  characteristic  would 
have  been  eliminated  by  low-pass  filtering  for  the  960-Hz  noise. 

31.  The  subroutine  containing  zero-  and  single-frequency  rejection 
filters  is  applicable  to  noise  frequencies  which  are  periodic  over  the 
time  history  to  be  plotted.  Its  usefulness  for  particular  jobs  should 

be  determined  from  examination  of  frequency  spectra  and  response  curves. 
Rejection  notches  may  not  be  overlapped  where  any  frequency  between  the 
intended  frequencies  must  be  retained  at  full  amplitude . 

Exercise  of  Judgment 

32.  Digital  filters  are  quite  powerful.  Nevertheless,  they  are  like 
other  tools  in  that  they  may  be  misused.  Though  the  filters  here  discussed 
are  applicable  to  digitized  time  histories  in  general,  their  blind  applica¬ 
tion  to  all  data  will  only  waste  time  and  money.  Professional  judgment 
must  guide  filter  use. 
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0.18  minute  per  thousand  samples 


EXAMPLE  SET  3,  FILTERED  DATA 

PLATE  8 
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TIME  FROH  DEI  -  SECS 


Gage  type:  Stress 

Sampling  frequency:  12,000  samples  per  second 

Filters  applied:  None 


EXAMPLE  SET  4,  RAW  DATA 


plait:  9 


impulse;  -  PS.I  X  SEC  STRESS  -  PS, I 

0.  5.  II).  IS.  -100  0  \  00  500 


SP3  SHI 

El  OCV 


f 


j 


l - 1 - —  ■» - 1 - 1 - r - - -i — 

-0.03  0.0’.  0.02  0.03  0.04  0.05  O.OS  0.07 

T THE  FROM  OET  -  SECS 
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Gage  type: 

Sampling  frequency: 
Filter  applied: 
Time  in  filter: 


Stress 

12,000  samples  per  second 
Single  frequency  rejection  (960hz) 
0.25  minute  per  thousand  samples 


EXAMPLE  SET  4,  FILTERED  DATA 


PLATE  .10 
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IS.  SRONSORtNO  MILITARY  ACTIVITY 


digital  tapes  produced  from  analog  field  data  often  carry  over  severe  noise  problems . 

A  batch  data  processing  capability  should  therefore  be  complemented  by  digital  filters 
which  can  remove  noise  routinely.  This  paper  describes  a  set  of  filters  developed  for 
use  with  Nuclear  Weapons  EffectB  Division  standard  data  processing  codes.  Emphasis  in 
assenfcly  of  the  set  was  placed  on  accuracy  of  signal  retention  and  on  adaptability  to 
general  purpose  -computers .  Part  I  discusses  data  which  is  sampled  at  equally  spaced 
times,  the  use  of  sample  numbers  to  designate  time  positions,  and  the  requirement  for 
use  of  a  sampling,  frequency  at  least  eight  times  greater  than  the  highest  expected 
signal  frequency.  Parts  II,  III,  and  IV  document  the  filters  of  the  three  subroutines 
presently  being  used  for  noise  removal.  The  first  subroutine  eliminates  randomly 
spaced  3ingle-sample  "spikes”  through  the  application  of  inequality  conditions  to 
data.  The  second  subroutine  provides  frequency  filters  for  the  removal  of  zero  drift 
and  single  noise  frequencies.  Lew- pass  frequency  filters  are  available  in  the  third 
subroutine.  Part  V  discusses  the  ten  example  plates  used  to  demonstrate  the  filter 
set's  effectiveness  end  speed  of  operation.  . 
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